Supplementary material D for ‘Incorporating unobserved heterogeneity in Weibull survival models: A Bayesian approach’
2017-03-06
Introduction
This vignette contains the code used for the case studies described in ‘Incorporating unobserved heterogeneity in Weibull survival models: A Bayesian approach’. All analyses were performed in R using the RMWreg library, which accompanies this manuscript.
Before starting the analysis, we proceed to open a new session of R Studio and to clean the existing R environment using
rm(list = ls())As a second step, we install and load the RMWreg.
library(RMWreg)
library(parallel) # For parallel run of multiple MCMC chains
library(coda) # To assess convergence of MCMC chains
library(gplots) # To us 'plotCI'
library(compiler); enableJIT(3) # For faster loops## [1] 0library(data.table)In the meantime, the RMWreg can be manually installed using the provided source file. However, the installation process will be simplified once this library is available through the CRAN repository.
Analysis of Autologous and Allogeneic Bone Marrow Transplant dataset
This dataset is publicly available as part of the KMsurv library. It contains post-surgery disease-free survival times (until relapse or death, in months) for 101 advanced acute myelogenous leukemia patients (51 right-censored observations). In the trial, 51 patients received an autologous bone marrow transplant, replacing the patient’s marrow with their own marrow treated with high doses of chemotherapy.The remaining patients had an allogeneic bone transplant, using marrow extracted from a sibling. Only the type of treatment is available as a covariate and thus an important amount of unobserved heterogeneity is expected.
Loading the data and creation of design matrix
library(KMsurv)
data(alloauto)
# Survival time and censoring indicator
Time = alloauto$time # Survival time
Event = alloauto$delta # 1: if observed event; 0: if censored
# Design matrix construction
x0 = rep(1, times = nrow(alloauto)) # Intercept
x1 = alloauto$type-1 # Treatment type - 0: allogeneic; 1: autologous
DesignMat = cbind(x0, x1) These data contain subjects.
Descriptive analysis
Number of subjects according to treatment type and censoring status
table(Event, x1)## x1
## Event 0 1
## 0 28 23
## 1 22 28Distribution of follow-up times (relapse, death or censoring) by treatment type
boxplot(Time ~ x1,
names = c("Allogeneic", "Autologous"),
xlab = "Treatment type", ylab = "Follow-up time",
frame = FALSE)
Distribution of observed event times (relapse or death, excluding censored observations) by treatment type
boxplot(Time[Event == 1] ~ x1[Event == 1],
names = c("Allogeneic", "Autologous"),
xlab = "Treatment type", ylab = "Observed event time",
frame = FALSE)
Fitting the model
Here, we analyse the data described above using Weibull model (no mixing) and RME-AFT regression models (i.e. RMW-AFT with \(\gamma = 0\)) using all the mixing distributions described in Table 1 of the main manuscript. The following function is used in order to run multiple chains in parallel.
RunRMWreg_MCMC <- function(Row, ParamTable,
Mixing, BaseModel,
N, Thin, Burn, Time, Event, DesignMat,
PrintProgress, StoreAdapt, lambdaPeriod,
StoreChains = FALSE, StoreDir = getwd(),
DataName)
{
require(RMWreg)
set.seed(ParamTable$seed[Row])
Chain <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat,
Mixing = Mixing, BaseModel = BaseModel,
PrintProgress = PrintProgress,
PriorCV = as.character(ParamTable$PriorCV[Row]),
PriorMeanCV = ParamTable$PriorMeanCV[Row],
StoreAdapt = StoreAdapt,
lambdaPeriod = lambdaPeriod,
StoreChains = StoreChains,
RunName = paste0(DataName,"_",Mixing,Row),
StoreDir = StoreDir)
return(Chain)
}Additionally, the following function is used to run graphical and formal convergence diagnostics for MCMC chains
PlotDiag <- function(Chain, titles)
{
for(i in 1:ncol(Chain))
{
plot(Chain[,i], main = titles[i],
cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5)
}
cat("Geweke diagnostics Z scores \n")
print(round(geweke.diag(Chain)$z,2)); cat("\n")
cat("Heidelberger and Welch's convergence diagnostic \n")
print(round(heidel.diag(Chain)[,2:3],2)); cat("\n")
cat("Effective Sample Size \n")
print(round(effectiveSize(Chain))); cat("\n")
}Throughout the following parameters are used
N = 600000; Thin = 50; Burn = 150000
StoreDir = "~/Documents/OneDrive/Projects/Survival/RMW/SuppMaterial/Chains"
titles = c(expression(beta[0]), expression(beta[1]), expression(gamma), expression(theta))Exponential (no mixing)
ParamTableNone = data.frame(seed = 1)set.seed(ParamTableNone$seed)
ChainNone <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat,
Mixing = "None", BaseModel = "Exponential",
PrintProgress = FALSE, StoreAdapt = TRUE,
StoreChains = TRUE, StoreDir = StoreDir,
RunName = "AA_None")ChainNone_mcmc = mcmc(ChainNone$beta)
colnames(ChainNone_mcmc) = titles[1:2]
PlotDiag(ChainNone_mcmc, titles = titles[1:2])
## Geweke diagnostics Z scores
## beta[0] beta[1]
## 0.26 -0.45
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.85
## beta[1] 1 0.65
##
## Effective Sample Size
## beta[0] beta[1]
## 9000 9000RME Exponential mixing
ParamTableExp = data.frame(seed = 2)set.seed(ParamTableExp$seed)
ChainExp <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat,
Mixing = "Exponential", BaseModel = "Exponential",
PrintProgress = FALSE, StoreAdapt = TRUE,
StoreChains = TRUE, StoreDir = StoreDir,
RunName = "AA_Exp", lambdaPeriod = 1)ChainExp_mcmc = mcmc(ChainExp$beta)
colnames(ChainExp_mcmc) = titles[1:2]
PlotDiag(ChainExp_mcmc, titles = titles[1:2])
## Geweke diagnostics Z scores
## beta[0] beta[1]
## -0.51 -0.25
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.37
## beta[1] 1 0.07
##
## Effective Sample Size
## beta[0] beta[1]
## 8062 9268RME Gamma mixing
ParamTableGam <- expand.grid("PriorCV" = c("TruncExp", "Pareto"),
"PriorMeanCV" = c(1.25, 1.5, 2, 5, 10))
set.seed(3)
ParamTableGam$seed = sample(1:1e6, size = nrow(ParamTableGam))
cat("Parameter values \n")## Parameter valuesprint(ParamTableGam); cat("\n")## PriorCV PriorMeanCV seed
## 1 TruncExp 1.25 168042
## 2 Pareto 1.25 807516
## 3 TruncExp 1.50 384942
## 4 Pareto 1.50 327734
## 5 TruncExp 2.00 602099
## 6 Pareto 2.00 604392
## 7 TruncExp 5.00 124633
## 8 Pareto 5.00 294599
## 9 TruncExp 10.00 577606
## 10 Pareto 10.00 630974cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainGam <- parLapply(cl, as.list(1:nrow(ParamTableGam)),
fun = RunRMWreg_MCMC,
ParamTable = ParamTableGam,
Mixing = "Gamma", BaseModel = "Exponential",
N = N, Thin = Thin, Burn = Burn,
Time = Time, Event = Event, DesignMat = DesignMat,
PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
StoreChains = TRUE, StoreDir = StoreDir, DataName = "AA")
stopCluster(cl)ChainGam_mcmc = lapply(ChainGam, function(x) {mcmc(cbind(x$beta, x$theta))})
ChainGam_mcmc = lapply(ChainGam_mcmc, function(x)
{
colnames(x) <- titles[c(1:2,4)]
return(x)
})
Plot = lapply(ChainGam_mcmc, PlotDiag, titles = titles[c(1:2,4)])
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 2.54 -1.64 0.57
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.08
## beta[1] 1 0.72
## theta 1 0.94
##
## Effective Sample Size
## beta[0] beta[1] theta
## 7296 8346 174
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -1.70 0.22 -0.88
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.30
## beta[1] 1 0.96
## theta 1 0.84
##
## Effective Sample Size
## beta[0] beta[1] theta
## 5625 9000 41
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.28 -0.77 -0.74
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.59
## beta[1] 1 0.62
## theta 1 0.90
##
## Effective Sample Size
## beta[0] beta[1] theta
## 5184 9000 39
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.27 0.22 -1.03
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.36
## beta[1] 1 0.59
## theta 1 0.82
##
## Effective Sample Size
## beta[0] beta[1] theta
## 7633 9000 918
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.15 0.78 -0.84
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.95
## beta[1] 1 0.29
## theta 1 0.71
##
## Effective Sample Size
## beta[0] beta[1] theta
## 5558 9000 115
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.53 0.79 1.34
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.39
## beta[1] 1 0.63
## theta 2701 0.42
##
## Effective Sample Size
## beta[0] beta[1] theta
## 8126 9000 1072
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.88 0.05 -1.20
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.08
## beta[1] 1 0.33
## theta 1 0.08
##
## Effective Sample Size
## beta[0] beta[1] theta
## 8643 8705 1550
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.19 0.35 -0.63
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.80
## beta[1] 1 0.51
## theta 1 0.99
##
## Effective Sample Size
## beta[0] beta[1] theta
## 8757 9000 2180
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 1.17 -0.42 0.53
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 3601 0.11
## beta[1] 1 0.49
## theta 1 0.13
##
## Effective Sample Size
## beta[0] beta[1] theta
## 4535 7861 89
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 0.14 0.52 -0.39
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.95
## beta[1] 1 0.38
## theta 1 0.82
##
## Effective Sample Size
## beta[0] beta[1] theta
## 7950 9000 1074RME Inverse Gamma mixing
ParamTableIGam <- expand.grid("PriorCV" = c("TruncExp", "Pareto"),
"PriorMeanCV" = c(1.25, 1.5))
set.seed(4)
ParamTableIGam$seed = sample(1:1e6, size = nrow(ParamTableIGam))
cat("Parameter values \n")## Parameter valuesprint(ParamTableIGam); cat("\n")## PriorCV PriorMeanCV seed
## 1 TruncExp 1.25 585801
## 2 Pareto 1.25 8946
## 3 TruncExp 1.50 293740
## 4 Pareto 1.50 277375cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainIGam <- parLapply(cl, as.list(1:nrow(ParamTableIGam)),
fun = RunRMWreg_MCMC,
ParamTable = ParamTableIGam,
Mixing = "InvGamma", BaseModel = "Exponential",
N = N, Thin = Thin, Burn = Burn,
Time = Time, Event = Event, DesignMat = DesignMat,
PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 5,
StoreChains = TRUE, StoreDir = StoreDir, DataName = "AA")
stopCluster(cl)ChainIGam_mcmc = lapply(ChainIGam, function(x) {mcmc(cbind(x$beta, x$theta))})
ChainIGam_mcmc = lapply(ChainIGam_mcmc, function(x)
{
colnames(x) <- titles[c(1:2,4)]
return(x)
})
Plot = lapply(ChainIGam_mcmc, PlotDiag, titles = titles[c(1:2,4)])
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 1.34 -1.02 -1.32
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.38
## beta[1] 1 0.34
## theta 1 0.53
##
## Effective Sample Size
## beta[0] beta[1] theta
## 1207 8238 609
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.42 -0.84 -0.29
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.80
## beta[1] 1 0.96
## theta 1 0.94
##
## Effective Sample Size
## beta[0] beta[1] theta
## 600 8264 180
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -1.42 0.19 0.99
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.26
## beta[1] 1 0.77
## theta 901 0.68
##
## Effective Sample Size
## beta[0] beta[1] theta
## 1371 8254 463
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 1.03 -0.38 -1.10
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.38
## beta[1] 1 0.82
## theta 1 0.22
##
## Effective Sample Size
## beta[0] beta[1] theta
## 1245 8592 418RME Inverse Gaussian mixing
ParamTableIGauss <- expand.grid("PriorCV" = c("TruncExp", "Pareto"),
"PriorMeanCV" = c(1.25, 1.5, 2))
set.seed(8)
ParamTableIGauss$seed = sample(1:1e6, size = nrow(ParamTableIGauss))
cat("Parameter values \n")## Parameter valuesprint(ParamTableIGauss); cat("\n")## PriorCV PriorMeanCV seed
## 1 TruncExp 1.25 466296
## 2 Pareto 1.25 207824
## 3 TruncExp 1.50 799657
## 4 Pareto 1.50 651870
## 5 TruncExp 2.00 321508
## 6 Pareto 2.00 718924cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainIGauss <- parLapply(cl, as.list(1:nrow(ParamTableIGauss)),
fun = RunRMWreg_MCMC,
ParamTable = ParamTableIGauss,
Mixing = "InvGauss", BaseModel = "Exponential",
N = N, Thin = Thin, Burn = Burn,
Time = Time, Event = Event, DesignMat = DesignMat,
PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 5,
StoreChains = TRUE, StoreDir = StoreDir, DataName = "AA")
stopCluster(cl)ChainIGauss_mcmc = lapply(ChainIGauss, function(x) {mcmc(cbind(x$beta, x$theta))})
ChainIGauss_mcmc = lapply(ChainIGauss_mcmc, function(x)
{
colnames(x) <- titles[c(1:2,4)]
return(x)
})
Plot = lapply(ChainIGauss_mcmc, PlotDiag, titles = titles[c(1:2,4)])
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.57 0.56 -0.73
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.69
## beta[1] 1 0.50
## theta 1 0.45
##
## Effective Sample Size
## beta[0] beta[1] theta
## 870 7752 1150
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 0.12 -0.60 -1.02
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.86
## beta[1] 1 0.66
## theta 1 0.91
##
## Effective Sample Size
## beta[0] beta[1] theta
## 448 7242 64
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 0.22 -0.38 -0.99
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.31
## beta[1] 1 0.99
## theta 1 0.14
##
## Effective Sample Size
## beta[0] beta[1] theta
## 1356 7032 391
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.58 -0.51 -1.24
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 3601 0.36
## beta[1] 1 0.80
## theta 1 0.05
##
## Effective Sample Size
## beta[0] beta[1] theta
## 394 7570 1273
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 0.61 0.74 0.65
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.84
## beta[1] 1 0.74
## theta 1 0.92
##
## Effective Sample Size
## beta[0] beta[1] theta
## 2062 7431 1886
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 0.94 -2.17 0.59
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.36
## beta[1] 1 0.30
## theta 1 0.64
##
## Effective Sample Size
## beta[0] beta[1] theta
## 1709 7579 1327RME Log-Normal mixing
ParamTableLN <- expand.grid("PriorCV" = c("TruncExp", "Pareto"),
"PriorMeanCV" = c(1.25, 1.5, 2, 5, 10))
set.seed(6)
ParamTableLN$seed = sample(1:1e6, size = nrow(ParamTableLN))
cat("Parameter values \n")## Parameter valuesprint(ParamTableLN); cat("\n")## PriorCV PriorMeanCV seed
## 1 TruncExp 1.25 606269
## 2 Pareto 1.25 937642
## 3 TruncExp 1.50 264352
## 4 Pareto 1.50 380093
## 5 TruncExp 2.00 807481
## 6 Pareto 2.00 978071
## 7 TruncExp 5.00 957928
## 8 Pareto 5.00 762727
## 9 TruncExp 10.00 509645
## 10 Pareto 10.00 64477cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainLN <- parLapply(cl, as.list(1:nrow(ParamTableLN)),
fun = RunRMWreg_MCMC,
ParamTable = ParamTableLN,
Mixing = "LogNormal", BaseModel = "Exponential",
N = N, Thin = Thin, Burn = Burn,
Time = Time, Event = Event, DesignMat = DesignMat,
PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 5,
StoreChains = TRUE, StoreDir = StoreDir, DataName = "AA")
stopCluster(cl)ChainLN_mcmc = lapply(ChainLN, function(x) {mcmc(cbind(x$beta, x$theta))})
ChainLN_mcmc = lapply(ChainLN_mcmc, function(x)
{
colnames(x) <- titles[c(1:2,4)]
return(x)
})
Plot = lapply(ChainLN_mcmc, PlotDiag, titles = titles[c(1:2,4)])
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.44 1.41 1.36
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.81
## beta[1] 1 0.09
## theta NA 0.02
##
## Effective Sample Size
## beta[0] beta[1] theta
## 7647 6243 946
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 1.62 -2.04 -0.58
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.54
## beta[1] 901 0.11
## theta 1 0.17
##
## Effective Sample Size
## beta[0] beta[1] theta
## 6217 5807 1169
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 0.12 -0.04 -0.94
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.78
## beta[1] 1 0.62
## theta 1 0.58
##
## Effective Sample Size
## beta[0] beta[1] theta
## 6860 6638 2093
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.12 -0.48 -0.02
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.47
## beta[1] 1 0.74
## theta 1 0.83
##
## Effective Sample Size
## beta[0] beta[1] theta
## 5942 5994 1463
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.91 1.36 -0.22
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.80
## beta[1] 1 0.78
## theta 1 0.27
##
## Effective Sample Size
## beta[0] beta[1] theta
## 5570 5478 2284
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.32 0.73 0.05
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.69
## beta[1] 1 0.13
## theta 1 0.55
##
## Effective Sample Size
## beta[0] beta[1] theta
## 5105 5001 1548
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.91 0.32 0.15
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.78
## beta[1] 1 0.59
## theta 1 0.14
##
## Effective Sample Size
## beta[0] beta[1] theta
## 4923 4847 2589
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## 0.99 -1.39 0.82
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.23
## beta[1] 1 0.74
## theta 1 0.53
##
## Effective Sample Size
## beta[0] beta[1] theta
## 4741 4810 1893
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.26 0.92 -0.56
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.29
## beta[1] 1 0.12
## theta 1 0.67
##
## Effective Sample Size
## beta[0] beta[1] theta
## 4317 4617 2802
## Geweke diagnostics Z scores
## beta[0] beta[1] theta
## -0.77 2.19 0.64
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.34
## beta[1] 1 0.29
## theta 1 0.16
##
## Effective Sample Size
## beta[0] beta[1] theta
## 4339 5005 1659Weibull (no mixing)
ParamTableWei = data.frame(seed = 7)set.seed(ParamTableWei$seed)
ChainWei <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat,
Mixing = "None", BaseModel = "Weibull",
Hyp1Gam = 4, Hyp2Gam = 1,
PrintProgress = FALSE, StoreAdapt = TRUE,
StoreChains = TRUE, StoreDir = StoreDir,
RunName = "AA_Wei")ChainWei_mcmc = mcmc(cbind(ChainWei$beta, ChainWei$gam))
colnames(ChainWei_mcmc) = titles[1:3]
PlotDiag(ChainWei_mcmc, titles = titles[1:3])
## Geweke diagnostics Z scores
## beta[0] beta[1] gamma
## -0.68 1.08 0.16
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.85
## beta[1] 1 0.92
## gamma 1 0.90
##
## Effective Sample Size
## beta[0] beta[1] gamma
## 9000 8660 9368RMW Exponential mixing
ParamTableRMWexp = data.frame(seed = 8)set.seed(ParamTableRMWexp$seed)
ChainRMWexp <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat,
Mixing = "Exponential", BaseModel = "Weibull",
Hyp1Gam = 4, Hyp2Gam = 1,
PrintProgress = FALSE, StoreAdapt = TRUE,
StoreChains = TRUE, StoreDir = StoreDir,
RunName = "AA_RMWexp", lambdaPeriod = 1)ChainRMWexp_mcmc = mcmc(cbind(ChainRMWexp$beta, ChainRMWexp$gam))
colnames(ChainRMWexp_mcmc) = titles[1:3]
PlotDiag(ChainRMWexp_mcmc, titles[1:3])
## Geweke diagnostics Z scores
## beta[0] beta[1] gamma
## 0.23 0.29 -0.23
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.31
## beta[1] 1 0.07
## gamma 1 0.77
##
## Effective Sample Size
## beta[0] beta[1] gamma
## 8472 8379 10191Model comparison
Pseudo ML
CDnone = RMWreg_CaseDeletion(ChainNone, Time, Event, DesignMat,
Mixing = "None", BaseModel = "Exponential")
ParamTableNone$logPsML <- sum(CDnone[,1])
CDexp = RMWreg_CaseDeletion(ChainExp, Time, Event, DesignMat,
Mixing = "Exponential", BaseModel = "Exponential")
ParamTableExp$logPsML <- sum(CDexp[,1])
CDgam = lapply(ChainGam, FUN = RMWreg_CaseDeletion,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "Gamma", BaseModel = "Exponential")
ParamTableGam$logPsML = sapply(CDgam, FUN = function(x) sum(x[,1]))
CDigam = lapply(ChainIGam, FUN = RMWreg_CaseDeletion,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "InvGamma", BaseModel = "Exponential")
ParamTableIGam$logPsML = sapply(CDigam, FUN = function(x) sum(x[,1]))
CDigauss = lapply(ChainIGauss, FUN = RMWreg_CaseDeletion,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "InvGauss", BaseModel = "Exponential")
ParamTableIGauss$logPsML = sapply(CDigauss, FUN = function(x) sum(x[,1]))
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
CDln = parLapply(cl, ChainLN, fun = RMWreg_CaseDeletion,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "LogNormal", BaseModel = "Exponential")
stopCluster(cl)
ParamTableLN$logPsML = sapply(CDln, FUN = function(x) sum(x[,1]))
CDwei = RMWreg_CaseDeletion(ChainWei, Time, Event, DesignMat,
Mixing = "None", BaseModel = "Weibull")
ParamTableWei$logPsML <- sum(CDwei[,1])
CDrmwexp = RMWreg_CaseDeletion(ChainRMWexp, Time, Event, DesignMat,
Mixing = "Exponential", BaseModel = "Weibull")
ParamTableRMWexp$logPsML <- sum(CDrmwexp[,1])ML
RunRMWreg_logML <- function(Row, ParamTable, Chains,
Mixing, BaseModel,
Time, Event, DesignMat,
Hyp1Gam = 1, Hyp2Gam = 1, N, Thin, Burn)
{
require(RMWreg)
logML = RMWreg_logML(Chains[[Row]], Time, Event, DesignMat,
PriorCV = as.character(ParamTable$PriorCV[Row]),
PriorMeanCV = ParamTable$PriorMeanCV[Row],
Hyp1Gam = Hyp1Gam, Hyp2Gam = Hyp2Gam,
Mixing = Mixing, BaseModel = BaseModel, N = N, Thin = Thin, Burn = Burn)
return(logML)
}
ParamTableNone$logML = RMWreg_logML(ChainNone, Time, Event, DesignMat,
Mixing = "None", BaseModel = "Exponential",
N = N, Thin = Thin, Burn = Burn)
ParamTableExp$logML = RMWreg_logML(ChainExp, Time, Event, DesignMat,
Mixing = "Exponential", BaseModel = "Exponential",
N = N, Thin = Thin, Burn = Burn)
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableGam$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableGam)),
fun = RunRMWreg_logML,
ParamTable = ParamTableGam,
Chains = ChainGam,
Mixing = "Gamma", BaseModel = "Exponential",
Time = Time, Event = Event, DesignMat = DesignMat,
N = N, Thin = Thin, Burn = Burn))
stopCluster(cl)
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableIGam$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableIGam)),
fun = RunRMWreg_logML,
ParamTable = ParamTableIGam,
Chains = ChainIGam,
Mixing = "InvGamma", BaseModel = "Exponential",
Time = Time, Event = Event, DesignMat = DesignMat,
N = N, Thin = Thin, Burn = Burn))
stopCluster(cl)
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableIGauss$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableIGauss)),
fun = RunRMWreg_logML,
ParamTable = ParamTableIGauss,
Chains = ChainIGauss,
Mixing = "InvGauss", BaseModel = "Exponential",
Time = Time, Event = Event, DesignMat = DesignMat,
N = N, Thin = Thin, Burn = Burn))
stopCluster(cl)
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableLN$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableLN)),
fun = RunRMWreg_logML,
ParamTable = ParamTableLN,
Chains = ChainLN,
Mixing = "LogNormal", BaseModel = "Exponential",
Time = Time, Event = Event, DesignMat = DesignMat,
N = N, Thin = Thin, Burn = Burn))
stopCluster(cl)
ParamTableWei$logML = RMWreg_logML(ChainWei, Time, Event, DesignMat,
Mixing = "None", BaseModel = "Weibull",
Hyp1Gam = 4, Hyp2Gam = 1,
N = N, Thin = Thin, Burn = Burn)
ParamTableRMWexp$logML = RMWreg_logML(ChainWei, Time, Event, DesignMat,
Mixing = "Exponential", BaseModel = "Weibull",
Hyp1Gam = 4, Hyp2Gam = 1,
N = N, Thin = Thin, Burn = Burn)DIC
ParamTableNone$DIC <- RMWreg_DIC(ChainNone, Time, Event, DesignMat,
Mixing = "None", BaseModel = "Exponential")
ParamTableExp$DIC <- RMWreg_DIC(ChainExp, Time, Event, DesignMat,
Mixing = "Exponential", BaseModel = "Exponential")
ParamTableGam$DIC = sapply(ChainGam, FUN = RMWreg_DIC,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "Gamma", BaseModel = "Exponential")
ParamTableIGam$DIC = sapply(ChainIGam, FUN = RMWreg_DIC,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "InvGamma", BaseModel = "Exponential")
ParamTableIGauss$DIC = sapply(ChainIGauss, FUN = RMWreg_DIC,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "InvGauss", BaseModel = "Exponential")
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableLN$DIC = simplify2array(parLapply(cl, ChainLN, fun = RMWreg_DIC,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "LogNormal", BaseModel = "Exponential"))
stopCluster(cl)
ParamTableWei$DIC <- RMWreg_DIC(ChainWei, Time, Event, DesignMat,
Mixing = "None", BaseModel = "Weibull")
ParamTableRMWexp$DIC <- RMWreg_DIC(ChainWei, Time, Event, DesignMat,
Mixing = "Exponential", BaseModel = "Weibull")Summary
cat("\n RME no mixing \n")##
## RME no mixingParamTableNone## seed logPsML logML DIC
## 1 1 -230.4599 -229.4623 459.9963cat("\n RME Exponential mixing \n")##
## RME Exponential mixingParamTableExp## seed logPsML logML DIC
## 1 2 -223.5881 -222.0672 446.5538cat("\n RME Gamma mixing \n")##
## RME Gamma mixingParamTableGam## PriorCV PriorMeanCV seed logPsML logML DIC
## 1 TruncExp 1.25 168042 -228.3560 -227.4952 455.9293
## 2 Pareto 1.25 807516 -228.2476 -227.0317 455.7122
## 3 TruncExp 1.50 384942 -227.6544 -226.7201 454.6764
## 4 Pareto 1.50 327734 -227.5672 -226.8410 454.4988
## 5 TruncExp 2.00 602099 -227.2252 -226.1494 453.8678
## 6 Pareto 2.00 604392 -227.2272 -226.4229 453.8926
## 7 TruncExp 5.00 124633 -226.6552 -225.3898 452.6433
## 8 Pareto 5.00 294599 -226.9003 -225.9342 453.2963
## 9 TruncExp 10.00 577606 -226.5228 -225.2106 452.3818
## 10 Pareto 10.00 630974 -226.8806 -225.8763 453.2428cat("\n RME Inverse Gamma mixing \n")##
## RME Inverse Gamma mixingParamTableIGam## PriorCV PriorMeanCV seed logPsML logML DIC
## 1 TruncExp 1.25 585801 -224.5198 -224.6049 448.8781
## 2 Pareto 1.25 8946 -224.5838 -224.7460 449.0500
## 3 TruncExp 1.50 293740 -224.2709 -223.8656 448.3744
## 4 Pareto 1.50 277375 -224.2931 -224.0833 448.4058cat("\n RME Inverse Gaussian mixing \n")##
## RME Inverse Gaussian mixingParamTableIGauss## PriorCV PriorMeanCV seed logPsML logML DIC
## 1 TruncExp 1.25 466296 -223.8578 -224.2714 447.7325
## 2 Pareto 1.25 207824 -223.9478 -224.2283 448.0506
## 3 TruncExp 1.50 799657 -223.4771 -223.2344 447.0858
## 4 Pareto 1.50 651870 -223.6137 -223.4765 447.3825
## 5 TruncExp 2.00 321508 -223.3694 -222.6429 446.8220
## 6 Pareto 2.00 718924 -223.3786 -222.9020 446.8793cat("\n RME Log-Normal mixing \n")##
## RME Log-Normal mixingParamTableLN## PriorCV PriorMeanCV seed logPsML logML DIC
## 1 TruncExp 1.25 606269 -225.2078 -225.5257 449.9074
## 2 Pareto 1.25 937642 -224.3059 -225.0513 448.2485
## 3 TruncExp 1.50 264352 -223.7999 -223.7006 447.3188
## 4 Pareto 1.50 380093 -223.4678 -223.5243 446.8311
## 5 TruncExp 2.00 807481 -223.0690 -222.3484 446.0267
## 6 Pareto 2.00 978071 -223.2085 -222.6695 446.3980
## 7 TruncExp 5.00 957928 -222.7849 -221.2237 445.6284
## 8 Pareto 5.00 762727 -223.0131 -222.0602 446.0679
## 9 TruncExp 10.00 509645 -222.9367 -221.2571 445.9855
## 10 Pareto 10.00 64477 -223.0865 -221.9649 446.2496cat("\n RMW no mixing \n")##
## RMW no mixingParamTableWei## seed logPsML logML DIC
## 1 7 -225.0809 -227.927 450.1554cat("\n RMW Exponential mixing \n")##
## RMW Exponential mixingParamTableRMWexp## seed logPsML logML DIC
## 1 8 -223.3197 -224.6084 450.69Figure 5
# Comparison Weibull vs Exponential
round(exp(ParamTableWei$logML - ParamTableNone$logML), 1)## [1] 4.6# Posterior summary for gam under Weibull model
round(median(ChainWei$gam),2)## [1] 0.69round(HPDinterval(mcmc(ChainWei$gam)),2)## lower upper
## V1 0.54 0.87
## attr(,"Probability")
## [1] 0.95# Comparison RMWexp vs RMEexp
round(exp(ParamTableExp$logML - ParamTableRMWexp$logML), 1)## [1] 12.7# Posterior summary for gam under Weibull model
round(median(ChainRMWexp$gam),2)## [1] 0.86round(HPDinterval(mcmc(ChainRMWexp$gam)),2)## lower upper
## V1 0.65 1.06
## attr(,"Probability")
## [1] 0.95# log-BFs
LOGBF1 = c(ParamTableNone$logML, ParamTableWei$logML, ParamTableExp$logML,
ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 1.25],
ParamTableIGam$logML[ParamTableIGam$PriorCV == "TruncExp" & ParamTableIGam$PriorMeanCV == 1.25],
ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 1.25],
ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 1.25]) -
ParamTableNone$logML
LOGBF1_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 1.25],
ParamTableIGam$logML[ParamTableIGam$PriorCV == "Pareto" & ParamTableIGam$PriorMeanCV == 1.25],
ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 1.25],
ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 1.25]) -
ParamTableNone$logML
LOGBF2 = c(ParamTableNone$logML, ParamTableWei$logML, ParamTableExp$logML,
ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 1.5],
ParamTableIGam$logML[ParamTableIGam$PriorCV == "TruncExp" & ParamTableIGam$PriorMeanCV == 1.5],
ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 1.5],
ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 1.5]) -
ParamTableNone$logML
LOGBF2_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 1.5],
ParamTableIGam$logML[ParamTableIGam$PriorCV == "Pareto" & ParamTableIGam$PriorMeanCV == 1.5],
ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 1.5],
ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 1.5]) -
ParamTableNone$logML
LOGBF3 = c(ParamTableNone$logML, ParamTableWei$logML, ParamTableExp$logML,
ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 2],
ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 2],
ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 2]) -
ParamTableNone$logML
LOGBF3_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 2],
ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 2],
ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 2]) -
ParamTableNone$logML
LOGBF4 = c(ParamTableNone$logML, ParamTableWei$logML, ParamTableExp$logML,
ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 5],
ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 5]) -
ParamTableNone$logML
LOGBF4_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 5],
ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 5]) -
ParamTableNone$logML
LOGBF5 = c(ParamTableNone$logML, ParamTableWei$logML, ParamTableExp$logML,
ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 10],
ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 10]) -
ParamTableNone$logML
LOGBF5_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 10],
ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 10]) -
ParamTableNone$logML
# log-PsBFs
LOGPsBF1 = c(ParamTableNone$logPsML, ParamTableWei$logPsML, ParamTableExp$logPsML,
ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 1.25],
ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "TruncExp" & ParamTableIGam$PriorMeanCV == 1.25],
ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 1.25],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 1.25]) -
ParamTableNone$logPsML
LOGPsBF1_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 1.25],
ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "Pareto" & ParamTableIGam$PriorMeanCV == 1.25],
ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 1.25],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 1.25]) -
ParamTableNone$logPsML
LOGPsBF2 = c(ParamTableNone$logPsML, ParamTableWei$logPsML, ParamTableExp$logPsML,
ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 1.5],
ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "TruncExp" & ParamTableIGam$PriorMeanCV == 1.5],
ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 1.5],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 1.5]) -
ParamTableNone$logPsML
LOGPsBF2_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 1.5],
ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "Pareto" & ParamTableIGam$PriorMeanCV == 1.5],
ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 1.5],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 1.5]) -
ParamTableNone$logPsML
LOGPsBF3 = c(ParamTableNone$logPsML, ParamTableWei$logPsML, ParamTableExp$logPsML,
ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 2],
ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 2],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 2]) -
ParamTableNone$logPsML
LOGPsBF3_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 2],
ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 2],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 2]) -
ParamTableNone$logPsML
LOGPsBF4 = c(ParamTableNone$logPsML, ParamTableWei$logPsML, ParamTableExp$logPsML,
ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 5],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 5]) -
ParamTableNone$logPsML
LOGPsBF4_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 5],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 5]) -
ParamTableNone$logPsML
LOGPsBF5 = c(ParamTableNone$logPsML, ParamTableWei$logPsML, ParamTableExp$logPsML,
ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 10],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 10]) -
ParamTableNone$logPsML
LOGPsBF5_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 10],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 10]) -
ParamTableNone$logPsML
par(mfrow=c(2,3))
par(mar=c(5,5,4,0.5))
plot(LOGBF1, LOGPsBF1, ylim = c(0,8.5), xlim = c(0,8.5),
ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
pch = c(4,3,8,0,5,1,2), cex = rep(3,7), bty = "n",
main = expression(paste("E(",c[v],")=1.25")),
cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd = 2)
points(LOGBF1_1, LOGPsBF1_1, pch = c(15,18,16,17), lwd = 2, cex = c(3,4.5,3,3))
plot(LOGBF2, LOGPsBF2, ylim = c(0,8.5), xlim = c(0,8.5),
ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
pch = c(4,3,8,0,5,1,2), cex = rep(3,7), bty = "n",
main = expression(paste("E(",c[v],")=1.5")),
cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd =2)
points(LOGBF2_1, LOGPsBF2_1, pch = c(15,18,16,17), lwd = 2, cex = c(3,4.5,3,3))
plot(LOGBF3, LOGPsBF3, ylim = c(0,8.5), xlim = c(0,8.5),
ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
pch = c(4,3,8,0,1,2), cex = rep(3,6), bty = "n",
main = expression(paste("E(",c[v],")=2")),
cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd = 2)
points(LOGBF3_1, LOGPsBF3_1, pch = c(15,16,17), lwd = 2, cex = rep(3,4))
plot(LOGBF4, LOGPsBF4, ylim = c(0,8.5), xlim = c(0,8.5),
ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
pch = c(4,3,8,0,2), cex = rep(3,5), bty = "n",
main = expression(paste("E(",c[v],")=5")),
cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd = 2)
points(LOGBF4_1, LOGPsBF4_1, pch = c(15,17), lwd = 2, cex = rep(3,4))
plot(LOGBF5, LOGPsBF5, ylim = c(0,8.5), xlim = c(0,8.5),
ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
pch = c(4,3,8,0,2), cex = rep(3,5), bty="n",
main = expression(paste("E(",c[v],")=10")),
cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd = 2)
points(LOGBF5_1, LOGPsBF5_1, pch = c(15,17), lwd = 2, cex = rep(3,4))
plot(1, type = "n", axes = F, xlab = "", ylab = "")
legend('center', c("Exponential","Weibull","RME-exp","RME-gam",
"RME-inv. gam","RME-inv. Gauss","RME-log normal"),
pch = c(4,3,8,0,5,1,2), col = "black", cex = 3, bty = "n")
Figure 6
prior = function(x, a, trunc, upper=Inf, type)
{
if(type=="TruncExp")
{
if(trunc==1){aux=dexp(x,rate=a)*exp(a)}
if(trunc==2){aux=dexp(x,rate=a)/(exp(-a)-exp(-upper*a))}
}
if(type=="Pareto")
{
if(trunc==1){aux=a*x^(-(a+1))}
if(trunc==2){aux=(a*x^(-(a+1)))/(1-upper^(-a))}
}
return(aux)
}
# CV as a function of theta
CV.RMEGAM<-function(theta) { sqrt(theta/(theta-2)) }
CV.RMEIGAM<-function(theta) { sqrt((theta+2)/(theta)) }
CV.RMEIGAUSS<-function(theta){ sqrt((5*theta^2+4*theta+1)/(theta^2+2*theta+1)) }
CV.RMELN<-function(theta){ sqrt(2*exp(theta)-1) }
par(mfrow=c(2,2))
par(mar=c(5, 5, 4, 2))
cv = seq(1, 50, by = 0.01)
hist(CV.RMELN(ChainLN[[7]]$theta), freq = FALSE,
ylim = c(0,0.4), xlim = c(1,10), breaks = 80,
main = expression(paste("Trunc. Exp. prior with ","E(",cv,")=5")),
ylab = "Frequency/Density", xlab = expression(R[cv]),
cex.lab = 2, cex.axis = 1.5, cex.main = 2.5)
lines(cv, prior(cv, a = 0.25, trunc = 1, upper = Inf, type = "TruncExp"), lwd = 4)
hist(CV.RMELN(ChainLN[[8]]$theta), freq = FALSE,
ylim = c(0,0.4), xlim = c(1,10), breaks = 100,
main = expression(paste("Pareto prior with ","E(",cv,")=5 ")),
ylab = "Frequency/Density", xlab = expression(R[cv]),
cex.lab = 2, cex.axis = 1.5, cex.main = 2.5)
lines(cv, prior(cv, a = 1.25, trunc = 1, upper = Inf, type = "Pareto"), lwd = 4)
hist(CV.RMELN(ChainLN[[9]]$theta), freq = FALSE,
ylim = c(0,0.3), xlim = c(1,10), breaks = 150,
main = expression(paste("Trunc. Exp. prior with ","E(",cv,")=10")),
ylab = "Frequency/Density", xlab = expression(R[cv]),
cex.lab = 2, cex.axis = 1.5, cex.main = 2.5)
lines(cv, prior(cv, a = 0.10, trunc = 1, upper = Inf, type = "TruncExp"), lwd = 4)
hist(CV.RMELN(ChainLN[[10]]$theta), freq = FALSE,
ylim = c(0,0.3), xlim = c(1,10), breaks = 200,
main = expression(paste("Pareto prior with ","E(",cv,")=10")),
ylab = "Frequency/Density", xlab = expression(R[cv]),
cex.lab = 2, cex.axis = 1.5, cex.main = 2.5)
lines(cv, prior(cv, a = 1.10, trunc = 1, upper = Inf, type = "Pareto"), lwd = 4)
Figure 7
This figure summarizes posterior estimates for the regression coefficients: posterior median and 95% Highest Posterior Density (HPD) interval
PM.b0 = c(median(ChainNone$beta[,1]),
median(ChainWei$beta[,1]),
median(ChainExp$beta[,1]),
unlist(lapply(ChainGam, function(x) median(x$beta[,1]))),
unlist(lapply(ChainIGam, function(x) median(x$beta[,1]))),
unlist(lapply(ChainIGauss, function(x) median(x$beta[,1]))),
unlist(lapply(ChainLN, function(x) median(x$beta[,1]))))
HPD.b0 = rbind(HPDinterval(ChainNone_mcmc[,1]),
HPDinterval(ChainWei_mcmc[,1]),
HPDinterval(ChainExp_mcmc[,1]),
t(sapply(ChainGam_mcmc, function(x) HPDinterval(x[,1]))),
t(sapply(ChainIGam_mcmc, function(x) HPDinterval(x[,1]))),
t(sapply(ChainIGauss_mcmc, function(x) HPDinterval(x[,1]))),
t(sapply(ChainLN_mcmc, function(x) HPDinterval(x[,1]))))
PM.b1 = c(median(ChainNone$beta[,2]),
median(ChainWei$beta[,2]),
median(ChainExp$beta[,2]),
unlist(lapply(ChainGam, function(x) median(x$beta[,2]))),
unlist(lapply(ChainIGam, function(x) median(x$beta[,2]))),
unlist(lapply(ChainIGauss, function(x) median(x$beta[,2]))),
unlist(lapply(ChainLN, function(x) median(x$beta[,2]))))
HPD.b1 = rbind(HPDinterval(ChainNone_mcmc[,2]),
HPDinterval(ChainWei_mcmc[,2]),
HPDinterval(ChainExp_mcmc[,2]),
t(sapply(ChainGam_mcmc, function(x) HPDinterval(x[,2]))),
t(sapply(ChainIGam_mcmc, function(x) HPDinterval(x[,2]))),
t(sapply(ChainIGauss_mcmc, function(x) HPDinterval(x[,2]))),
t(sapply(ChainLN_mcmc, function(x) HPDinterval(x[,2]))))
par(mfrow = c(2, 1))
par(mar = c(5,6,4,0.5))
plotCI(x = PM.b0, uiw = HPD.b0[,2] - PM.b0, liw = PM.b0 - HPD.b0[,1],
main = expression("Intercept"), xaxt = "n", lwd = 2,
xlab = "Mixing-Prior", ylab = expression(beta[0]),
cex.lab = 2, cex.axis = 2, cex.main = 2.5, gap = 1.5)
axis(side = 1, at = c(1.5,3,6,11,14.3,16.6,19,22,26,31),
labels = c("None", "Exp", "Gamma-Trunc.Exp.", "Gam-Pareto",
"IGam-T.E.", "IGam-P", "IGauss-T.E.", "I.Gauss-P.",
"LN-Trunc.Exp", "LN-Pareto"), cex.axis = 1)
abline(v = c(2.5,3.5,8.5,13.5,15.5,17.5,20.5,23.5,28.5), lty = 2)
abline(h = 0, lty = 3)
plotCI(x = PM.b1, uiw = HPD.b1[,2] - PM.b1, liw = PM.b1 - HPD.b1[,1],
main = expression("Treatment: autologous"), xaxt = "n", lwd = 2,
xlab = "Mixing-Prior", ylab = expression(beta[1]),
cex.lab = 2, cex.axis = 2, cex.main = 2.5, gap = 1.5)
axis(side = 1, at = c(1.5,3,6,11,14.3,16.6,19,22,26,31),
labels = c("None", "Exp", "Gamma-Trunc.Exp.", "Gam-Pareto",
"IGam-T.E.", "IGam-P", "IGauss-T.E.", "I.Gauss-P.",
"LN-Trunc.Exp", "LN-Pareto"), cex.axis = 1)
abline(v = c(2.5,3.5,8.5,13.5,15.5,17.5,20.5,23.5,28.5), lty = 2)
abline(h = 0, lty = 3)
Figure 8
BoxplotLambda <- function(Chain, main, ref, ylim, cex.lab, cex.axis, cex.main)
{
PM = as.vector(apply(Chain$lambda, 2, "median"))
U = as.vector(HPDinterval(mcmc(Chain$lambda))[,2]) - PM
L = PM - as.vector(HPDinterval(mcmc(Chain$lambda))[,1])
plotCI(x = PM, uiw = U, liw = L,
main = main, xaxt = "n", lwd = 1.5,
xlab = "Patient", ylab = expression(lambda[i]),
cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, ylim = ylim, cex = 1)
abline(h = ref, lty = 2)
return(PM)
}
# Outlier detection
BFexp = RMWreg_BFoutlier(ChainExp, RefLambda = Event + (1-Event)*1/2,
Time, Event, DesignMat,
Mixing = "Exponential", BaseModel = "Exponential")
par(mar = c(5,6,4,0.5))
par(mfrow = c(2,1))
PM.RME.EXP = BoxplotLambda(Chain = ChainExp, main = "(a)",
ref = c(0.5,1), ylim = c(0,5), cex.lab = 2, cex.axis = 2, cex.main=2)## Note: no visible binding for global variable 'lambda'
## Note: no visible binding for global variable 'i'points((1:length(Time))[Event == 0], PM.RME.EXP[Event == 0], lwd = 2, pch = 16)
axis(side = 1, at = c(25,75), labels = c("Treatment: Allogeneic","Treatment: Autologous"), cex.axis = 2)
plot(BFexp, type = "l", ylim = c(0,2), main = "(b)",
cex.main = 2, cex.lab = 2, cex.axis = 2, ylab = "Bayes Factors", xlab = "Patient")
abline(h = 1, lty = 2)
Analysis of Cerebral Palsy dataset
These data contain records of 1,549 children affected by cerebral palsy and born during the period 1966-1984 in the area of the Mersey Region Health Authority. The time to follow-up (survival or censoring) is recorded in years since birth. Following the analysis in Kwong and Hutton (2003), we use the amount of severe impairments (ambulation, manual dexterity and mental ability) and the birth weight (in kg.) as predictors for the time to death.
These data were kindly provided by P.O.D. Pharoah and Jane Hutton and we are not authorized to release it. However, this file contains relevant summary information and the code used to perform the analysis. These should serve as a guidance for potential users of the RMWreg library.
Loading the data and creation of design matrix
In this example, the data file has been stored in a data.path folder.
Data = read.csv(file.path(data.path, "CPdata.csv"), sep="")
# Removal of subjects with missing birth weight
Data = Data[is.na(Data$bw)==FALSE,]
# Survival time and censoring indicator
Time = Data$Lifeyr # Survival time
Event = Data$Dead # 1: if observed event; 0: if censored
# Design matrix construction
x0 = rep(1, times = nrow(Data)) # Intercept
x1 = Data$Ambul + Data$Mandex + Data$Mentab # No of impairments (treated as categorical)
x1.0 = as.numeric(I(x1 == 0))
x1.1 = as.numeric(I(x1 == 1))
x1.2 = as.numeric(I(x1 == 2))
x2 = Data$bw / 1000 # Birth weights in kgs
DesignMat = cbind(x0, x1.0, x1.1, x1.2, x2)These data contain 1549 subjects.
Descriptive analysis
Number of subjects according to the number of impairments and censoring status
table(Event, x1)## x1
## Event 0 1 2 3
## 0 960 208 49 90
## 1 23 25 21 173Number of subjects according to birth weight and censoring status
bw_group = ifelse(Data$bw < 1500, "1. Very low",
ifelse(Data$bw < 2500, "2. Low", "3. Normal"))
table(Event, bw_group)## bw_group
## Event 1. Very low 2. Low 3. Normal
## 0 197 341 769
## 1 24 56 162Distribution of follow-up times (survival or censoring) by the number of impairments
boxplot(Time ~ x1,
names = 0:3,
xlab = "No. of impairments", ylab = "Follow-up time",
frame = FALSE)
Distribution of observed event times (relapse or death, excluding censored observations) by the number of impairments
boxplot(Time[Event == 1] ~ x1[Event == 1],
names = 0:3,
xlab = "No. of impairments", ylab = "Observed event time",
frame = FALSE)
Distribution of follow-up times (survival or censoring) by birth weight group
boxplot(Time ~ bw_group,
xlab = "Birth weight group", ylab = "Follow-up time",
frame = FALSE)
Distribution of observed event times (relapse or death, excluding censored observations) by birth weight group
boxplot(Time[Event == 1] ~ bw_group[Event == 1],
xlab = "Birth weight group", ylab = "Observed event time",
frame = FALSE)
RMW-AFT analysis
RunRMWreg_MCMC <- function(Row, ParamTable,
Mixing, BaseModel,
N, Thin, Burn, Time, Event, DesignMat,
PrintProgress, StoreAdapt, lambdaPeriod,
StoreChains = FALSE, StoreDir = getwd(),
DataName)
{
require(RMWreg)
if("PriorCV" %in% names(ParamTable))
{
set.seed(ParamTable$seed[Row])
Chain <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat,
Mixing = Mixing, BaseModel = BaseModel,
PrintProgress = PrintProgress,
PriorCV = as.character(ParamTable$PriorCV[Row]),
PriorMeanCV = ParamTable$PriorMeanCV[Row],
Hyp1Gam = ParamTable$Hyp1Gam[Row],
Hyp2Gam = ParamTable$Hyp2Gam[Row],
StoreAdapt = StoreAdapt,
lambdaPeriod = lambdaPeriod,
StoreChains = StoreChains,
RunName = paste0(DataName,"_",Mixing,Row),
StoreDir = StoreDir)
}
else
{
set.seed(ParamTable$seed[Row])
Chain <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat,
Mixing = Mixing, BaseModel = BaseModel,
PrintProgress = PrintProgress,
Hyp1Gam = ParamTable$Hyp1Gam[Row],
Hyp2Gam = ParamTable$Hyp2Gam[Row],
StoreAdapt = StoreAdapt,
lambdaPeriod = lambdaPeriod,
StoreChains = StoreChains,
RunName = paste0(DataName,"_",Mixing,Row),
StoreDir = StoreDir)
}
return(Chain)
}Additionally, the following function is used to run graphical and formal convergence diagnostics for MCMC chains
PlotDiag <- function(Chain, titles)
{
for(i in 1:ncol(Chain))
{
plot(Chain[,i], main = titles[i],
cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5)
}
cat("Geweke diagnostics Z scores \n")
print(round(geweke.diag(Chain)$z,2)); cat("\n")
cat("Heidelberger and Welch's convergence diagnostic \n")
print(round(heidel.diag(Chain)[,2:3],2)); cat("\n")
cat("Effective Sample Size \n")
print(round(effectiveSize(Chain))); cat("\n")
}Throughout the following parameters are used
N = 600000; Thin = 50; Burn = 150000
StoreDir = "~/Documents/OneDrive/Projects/Survival/RMW/SuppMaterial/Chains"
titles = c(expression(beta[0]), expression(beta[1]),
expression(beta[2]), expression(beta[3]), expression(beta[4]),
expression(gamma), expression(theta))Weibull (no mixing)
ParamTableNone <- cbind.data.frame("Hyp1Gam" = c(4, 1), "Hyp2Gam" = c(1, 1))
set.seed(100)
ParamTableNone$seed = sample(1:1e6, size = nrow(ParamTableNone))
cat("Parameter values \n")## Parameter valuesprint(ParamTableNone); cat("\n")## Hyp1Gam Hyp2Gam seed
## 1 4 1 307767
## 2 1 1 257673cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainNone <- parLapply(cl, as.list(1:nrow(ParamTableNone)),
fun = RunRMWreg_MCMC,
ParamTable = ParamTableNone,
Mixing = "None", BaseModel = "Weibull",
N = N, Thin = Thin, Burn = Burn,
Time = Time, Event = Event, DesignMat = DesignMat,
PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)ChainNone_mcmc = lapply(ChainNone, function(x) {mcmc(cbind(x$beta, x$gam))})
ChainNone_mcmc = lapply(ChainNone_mcmc, function(x)
{
colnames(x) <- titles[1:6]
return(x)
})
Plot = lapply(ChainNone_mcmc, PlotDiag, titles = titles[1:6])
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma
## 1.33 0.19 0.43 1.03 -1.37 -1.01
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.66
## beta[1] 1 0.31
## beta[2] 1 0.94
## beta[3] 1 0.67
## beta[4] 1 0.44
## gamma 1 0.47
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma
## 3653 9000 9490 8667 3761 8594
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma
## -1.06 -0.71 0.74 -0.18 1.02 -0.05
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.55
## beta[1] 1 0.62
## beta[2] 1 0.50
## beta[3] 1 0.52
## beta[4] 1 0.59
## gamma 1 0.26
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma
## 3975 9000 9000 9000 3997 9000RMW Exponential mixing
ParamTableExp <- cbind.data.frame("Hyp1Gam" = c(4, 1), "Hyp2Gam" = c(1, 1))
set.seed(101)
ParamTableExp$seed = sample(1:1e6, size = nrow(ParamTableExp))
cat("Parameter values \n")## Parameter valuesprint(ParamTableExp); cat("\n")## Hyp1Gam Hyp2Gam seed
## 1 4 1 372199
## 2 1 1 43825cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainExp <- parLapply(cl, as.list(1:nrow(ParamTableExp)),
fun = RunRMWreg_MCMC,
ParamTable = ParamTableExp,
Mixing = "Exponential", BaseModel = "Weibull",
N = N, Thin = Thin, Burn = Burn,
Time = Time, Event = Event, DesignMat = DesignMat,
PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)ChainExp_mcmc = lapply(ChainExp, function(x) {mcmc(cbind(x$beta, x$gam))})
ChainExp_mcmc = lapply(ChainExp_mcmc, function(x)
{
colnames(x) <- titles[1:6]
return(x)
})
Plot = lapply(ChainExp_mcmc, PlotDiag, titles = titles[1:6])
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma
## 1.19 -2.67 -0.26 0.25 -0.98 1.17
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.36
## beta[1] 901 0.22
## beta[2] 1 0.76
## beta[3] 1 0.89
## beta[4] 1 0.42
## gamma 1 0.12
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma
## 2514 9000 9000 9000 2644 9000
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma
## 0.21 1.55 1.03 -1.82 -0.04 -1.13
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.67
## beta[1] 1 0.49
## beta[2] 1 0.27
## beta[3] 901 0.23
## beta[4] 1 0.61
## gamma 1 0.66
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma
## 2519 8387 9000 9000 2593 8545RMW Gamma mixing
ParamTableGam <- expand.grid("Hyp1Gam" = c(4, 1),
"PriorMeanCV" = c(1.5, 5),
"PriorCV" = c("TruncExp", "Pareto"))
ParamTableGam$Hyp2Gam = with(ParamTableGam,
ifelse(Hyp1Gam == 4, 1, ifelse(Hyp1Gam == 1, 1, 0.1)))
set.seed(102)
ParamTableGam$seed = sample(1:1e6, size = nrow(ParamTableGam))
cat("Parameter values \n")## Parameter valuesprint(ParamTableGam); cat("\n")## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed
## 1 4 1.5 TruncExp 1 571629
## 2 1 1.5 TruncExp 1 492883
## 3 4 5.0 TruncExp 1 783694
## 4 1 5.0 TruncExp 1 540340
## 5 4 1.5 Pareto 1 88002
## 6 1 1.5 Pareto 1 420628
## 7 4 5.0 Pareto 1 976328
## 8 1 5.0 Pareto 1 330313cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainGam <- parLapply(cl, as.list(1:nrow(ParamTableGam)),
fun = RunRMWreg_MCMC,
ParamTable = ParamTableGam,
Mixing = "Gamma", BaseModel = "Weibull",
N = 2*N, Thin = 2*Thin, Burn = 2*Burn,
Time = Time, Event = Event, DesignMat = DesignMat,
PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)ChainGam_mcmc = lapply(ChainGam, function(x) {mcmc(cbind(x$beta, x$gam, x$theta))})
ChainGam_mcmc = lapply(ChainGam_mcmc, function(x)
{
colnames(x) <- titles
return(x)
})
Plot = lapply(ChainGam_mcmc, PlotDiag, titles = titles)
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 0.04 0.37 -1.09 0.62 -0.12 0.91 -0.41
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.80
## beta[1] 1 0.65
## beta[2] 1 0.69
## beta[3] 1 0.99
## beta[4] 1 0.66
## gamma 1 0.22
## theta 1 0.97
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 5231 7794 8025 9000 5446 4483 1241
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 2.38 2.33 1.00 -0.17 -1.37 -1.71 1.43
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.18
## beta[1] 1 0.75
## beta[2] 1 0.62
## beta[3] 1 0.13
## beta[4] 1 0.36
## gamma 1 0.75
## theta 1 0.52
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 5110 7449 9000 9000 5453 3989 1248
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 0.97 1.36 0.90 -0.56 -0.72 -1.74 1.22
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.72
## beta[1] 1 0.27
## beta[2] 1 0.25
## beta[3] 1 0.38
## beta[4] 1 0.42
## gamma 1 0.16
## theta 1 0.20
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 4892 6488 8290 8695 5131 4037 2767
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -0.56 -0.60 0.56 0.77 0.23 0.56 0.58
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.18
## beta[1] 1 0.45
## beta[2] 1 0.76
## beta[3] 1 0.09
## beta[4] 1 0.28
## gamma 1 0.37
## theta 1 0.16
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 5138 6517 8150 9000 5225 4202 1330
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -1.55 0.48 1.56 -0.27 0.28 1.58 -1.79
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.58
## beta[1] 1 0.61
## beta[2] 1 0.24
## beta[3] 1 0.83
## beta[4] 1 0.85
## gamma 1 0.97
## theta 1 0.64
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 5213 7603 9000 9000 5615 3104 472
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -2.02 -1.62 -1.44 -1.73 0.22 3.80 -3.35
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.16
## beta[1] 1 0.59
## beta[2] 1 0.62
## beta[3] 1 0.30
## beta[4] 1 0.95
## gamma 901 0.38
## theta 1 0.60
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 4224 7568 9000 9000 5564 2169 204
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -0.53 -2.22 0.87 -0.70 0.18 1.16 -0.70
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.60
## beta[1] 1 0.22
## beta[2] 1 0.25
## beta[3] 1 0.49
## beta[4] 1 0.44
## gamma 1 0.83
## theta 1 0.54
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 5107 7594 9000 9000 5392 4375 1319
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -0.28 0.45 -1.01 0.56 0.70 -0.19 -0.01
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.37
## beta[1] 1 0.26
## beta[2] 1 0.80
## beta[3] 1 0.83
## beta[4] 1 0.53
## gamma 1 0.45
## theta 1 0.72
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 5024 7968 8951 9794 5237 4396 978RMW Inverse Gamma mixing
ParamTableIGam <- expand.grid("Hyp1Gam" = c(4, 1),
"PriorMeanCV" = c(1.5, 5),
"PriorCV" = c("TruncExp", "Pareto"))
ParamTableIGam$Hyp2Gam = with(ParamTableIGam,
ifelse(Hyp1Gam == 4, 1, ifelse(Hyp1Gam == 1, 1, 0.1)))
set.seed(103)
ParamTableIGam$seed = sample(1:1e6, size = nrow(ParamTableIGam))
cat("Parameter values \n")## Parameter valuesprint(ParamTableIGam); cat("\n")## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed
## 1 4 1.5 TruncExp 1 215942
## 2 1 1.5 TruncExp 1 63169
## 3 4 5.0 TruncExp 1 521826
## 4 1 5.0 TruncExp 1 503616
## 5 4 1.5 Pareto 1 120486
## 6 1 1.5 Pareto 1 87275
## 7 4 5.0 Pareto 1 433560
## 8 1 5.0 Pareto 1 193080cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainIGam <- parLapply(cl, as.list(1:nrow(ParamTableIGam)),
fun = RunRMWreg_MCMC,
ParamTable = ParamTableIGam,
Mixing = "InvGamma", BaseModel = "Weibull",
N = N * 2, Thin = Thin * 2, Burn = Burn * 2,
Time = Time, Event = Event, DesignMat = DesignMat,
PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)ChainIGam_mcmc = lapply(ChainIGam, function(x) {mcmc(cbind(x$beta, x$gam, x$theta))})
ChainIGam_mcmc = lapply(ChainIGam_mcmc, function(x)
{
colnames(x) <- titles
return(x)
})
Plot = lapply(ChainIGam_mcmc, PlotDiag, titles = titles)
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 1.16 -1.21 -1.46 -0.55 0.52 0.93 -1.20
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.27
## beta[1] 1 0.48
## beta[2] 1 0.21
## beta[3] 1 0.88
## beta[4] 1 0.60
## gamma 1 0.35
## theta 1 0.45
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 315 9000 7819 8420 4446 456 114
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 0.36 -0.79 -0.60 -0.46 -0.69 0.28 -0.29
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.91
## beta[1] 1 0.42
## beta[2] 1 0.19
## beta[3] 1 0.86
## beta[4] 1 0.32
## gamma 1 0.87
## theta 1 0.93
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 331 8629 9000 9000 4967 529 187
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -1.51 2.10 1.39 -0.13 0.58 -1.68 1.11
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1801 0.56
## beta[1] 901 0.19
## beta[2] 1 0.34
## beta[3] 1 0.43
## beta[4] 1 0.50
## gamma 1801 0.18
## theta 1801 0.73
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 360 9000 9000 9000 4849 427 193
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 0.93 -1.70 -1.02 1.53 0.91 1.05 -1.19
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.28
## beta[1] 1 0.23
## beta[2] 1 0.66
## beta[3] 1 0.19
## beta[4] 1 0.68
## gamma 1 0.44
## theta 1 0.20
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 241 8055 9000 8020 4962 461 146
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -0.74 0.13 -0.55 -1.01 -1.39 -0.80 0.77
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.25
## beta[1] 1 0.61
## beta[2] 1 0.60
## beta[3] 1 0.39
## beta[4] 1 0.61
## gamma 1 0.36
## theta 1 0.09
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 230 9000 9000 9000 4620 474 111
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -1.20 1.30 1.38 -0.35 0.21 -1.45 0.81
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.85
## beta[1] 1 0.49
## beta[2] 1 0.44
## beta[3] 1 0.66
## beta[4] 1 0.20
## gamma 1 0.48
## theta 1 0.86
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 104 5578 9000 9000 5000 267 63
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 1.64 0.44 0.49 0.01 -0.72 1.06 -1.16
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] NA 0.00
## beta[1] NA 0.00
## beta[2] 1 0.83
## beta[3] 1 0.84
## beta[4] 1 0.21
## gamma NA 0.00
## theta NA 0.00
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 42 5087 9000 8630 4750 251 12
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 0.96 -1.37 -1.31 -0.51 -0.43 0.94 -1.08
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1801 0.08
## beta[1] 1 0.52
## beta[2] 1 0.29
## beta[3] 1 0.74
## beta[4] 1 0.92
## gamma 1801 0.06
## theta 1 0.08
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 295 9000 9000 8105 5115 553 181RMW Inverse Gaussian mixing
ParamTableIGauss <- expand.grid("Hyp1Gam" = c(4, 1),
"PriorMeanCV" = c(1.5, 5),
"PriorCV" = c("TruncExp", "Pareto"))
ParamTableIGauss$Hyp2Gam = with(ParamTableIGauss,
ifelse(Hyp1Gam == 4, 1, ifelse(Hyp1Gam == 1, 1, 0.1)))
set.seed(104)
ParamTableIGauss$seed = sample(1:1e6, size = nrow(ParamTableIGauss))
cat("Parameter values \n")## Parameter valuesprint(ParamTableIGauss); cat("\n")## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed
## 1 4 1.5 TruncExp 1 364452
## 2 1 1.5 TruncExp 1 772498
## 3 4 5.0 TruncExp 1 734878
## 4 1 5.0 TruncExp 1 972528
## 5 4 1.5 Pareto 1 740140
## 6 1 1.5 Pareto 1 200713
## 7 4 5.0 Pareto 1 377307
## 8 1 5.0 Pareto 1 247017cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainIGauss <- parLapply(cl, as.list(1:nrow(ParamTableIGauss)),
fun = RunRMWreg_MCMC,
ParamTable = ParamTableIGauss,
Mixing = "InvGauss", BaseModel = "Weibull",
N = 2*N, Thin = 2*Thin, Burn = 2*Burn,
Time = Time, Event = Event, DesignMat = DesignMat,
PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)ChainIGauss_mcmc = lapply(ChainIGauss, function(x) {mcmc(cbind(x$beta, x$gam, x$theta))})
ChainIGauss_mcmc = lapply(ChainIGauss_mcmc, function(x)
{
colnames(x) <- titles
return(x)
})
Plot = lapply(ChainIGauss_mcmc, PlotDiag, titles = titles)
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 0.83 0.79 0.80 -0.77 -0.96 0.11 0.61
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.24
## beta[1] 1 0.07
## beta[2] 1 0.88
## beta[3] 1 0.95
## beta[4] 1 0.53
## gamma 1 0.19
## theta 1 0.06
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 298 5342 7007 7754 1096 217 134
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -0.70 0.25 0.52 1.06 0.36 -1.16 -1.60
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] NA 0.00
## beta[1] 1 0.09
## beta[2] 1 0.12
## beta[3] 1 0.66
## beta[4] 1 0.12
## gamma 1 0.06
## theta 1 0.15
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 296 3266 7750 8240 1040 142 115
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -0.03 0.38 -0.91 1.63 -0.18 0.26 0.05
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.05
## beta[1] 1 0.14
## beta[2] 1 0.91
## beta[3] 1 0.14
## beta[4] 1 0.79
## gamma 1 0.12
## theta 1801 0.07
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 843 6211 8506 9000 4218 531 229
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -1.06 1.23 0.37 -1.27 -0.53 -0.99 -0.69
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.27
## beta[1] 1 0.33
## beta[2] 1 0.61
## beta[3] 1 0.44
## beta[4] 1 0.51
## gamma 1 0.49
## theta 1 0.72
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 597 6137 8700 9000 3971 442 341
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 1.49 -0.86 -0.31 1.55 -0.78 0.42 -0.36
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.37
## beta[1] 1 0.36
## beta[2] 1 0.07
## beta[3] 1 0.46
## beta[4] 1 0.64
## gamma 1 0.49
## theta 1 0.69
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 478 5766 9000 9000 4171 469 424
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -0.51 0.32 -0.71 -0.13 0.31 -0.41 -0.26
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.90
## beta[1] 1 0.57
## beta[2] 1 0.20
## beta[3] 1 0.94
## beta[4] 1 0.52
## gamma 1 0.82
## theta 1 0.80
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 323 4792 9000 8370 4278 372 275
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 0.50 -0.12 -0.29 0.96 -1.17 0.14 -0.61
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.63
## beta[1] 1 0.48
## beta[2] 1 0.34
## beta[3] 1 0.32
## beta[4] 1 0.22
## gamma 1 0.49
## theta 1 0.63
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 559 6760 8362 9000 4113 447 360
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 2.12 -1.77 0.65 -0.23 0.81 2.17 1.33
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.26
## beta[1] 1 0.63
## beta[2] 1 0.40
## beta[3] 1 0.85
## beta[4] 1 0.60
## gamma 901 0.09
## theta 2701 0.44
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 372 5317 8534 9000 4180 310 194RMW Log Normal mixing
ParamTableLN <- expand.grid("Hyp1Gam" = c(4, 1),
"PriorMeanCV" = c(1.5, 5),
"PriorCV" = c("TruncExp", "Pareto"))
ParamTableLN$Hyp2Gam = with(ParamTableLN,
ifelse(Hyp1Gam == 4, 1, ifelse(Hyp1Gam == 1, 1, 0.1)))
cat("Parameter values \n")## Parameter valuesprint(ParamTableLN); cat("\n")## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam
## 1 4 1.5 TruncExp 1
## 2 1 1.5 TruncExp 1
## 3 4 5.0 TruncExp 1
## 4 1 5.0 TruncExp 1
## 5 4 1.5 Pareto 1
## 6 1 1.5 Pareto 1
## 7 4 5.0 Pareto 1
## 8 1 5.0 Pareto 1set.seed(105)
ParamTableLN$seed = sample(1:1e6, size = nrow(ParamTableLN))cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainLN <- parLapply(cl, as.list(1:nrow(ParamTableLN)),
fun = RunRMWreg_MCMC,
ParamTable = ParamTableLN,
Mixing = "LogNormal", BaseModel = "Weibull",
N = 2*N, Thin = 2*Thin, Burn = 2*Burn,
Time = Time, Event = Event, DesignMat = DesignMat,
PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)ChainLN_mcmc = lapply(ChainLN, function(x) {mcmc(cbind(x$beta, x$gam, x$theta))})
ChainLN_mcmc = lapply(ChainLN_mcmc, function(x)
{
colnames(x) <- titles
return(x)
})
Plot = lapply(ChainLN_mcmc, PlotDiag, titles = titles)
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 0.40 -1.89 -0.99 -0.84 -0.92 1.88 1.02
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.68
## beta[1] 1 0.07
## beta[2] 1 0.58
## beta[3] 1 0.07
## beta[4] 1 0.64
## gamma 1 0.13
## theta 1 0.24
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 3359 5908 8689 9000 3504 1044 980
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 0.16 0.38 0.13 -1.55 0.11 -0.36 0.01
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.82
## beta[1] 1 0.28
## beta[2] 1 0.25
## beta[3] 1 0.29
## beta[4] 1 0.88
## gamma 1 0.18
## theta 1 0.29
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 3418 4761 8380 9000 3581 841 750
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -0.87 -0.83 1.16 1.34 0.67 0.02 -0.16
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.37
## beta[1] 1 0.33
## beta[2] 1 0.79
## beta[3] 1 0.43
## beta[4] 1 0.36
## gamma 1 0.94
## theta 1 0.89
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 3157 6741 8267 9000 3302 1055 978
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 1.54 0.16 -0.71 0.82 -1.65 0.53 0.68
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.85
## beta[1] 1 0.13
## beta[2] 1 0.40
## beta[3] 1 0.19
## beta[4] 1 0.48
## gamma 1 0.49
## theta 1 0.69
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 2918 5410 8573 9000 3331 844 856
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -1.10 1.68 -0.59 -0.87 1.85 -2.09 -2.14
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.27
## beta[1] 1 0.50
## beta[2] 1 0.61
## beta[3] 1 0.70
## beta[4] 1 0.11
## gamma 1 0.13
## theta 1 0.07
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 3776 6280 8587 9000 3848 851 846
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -0.34 0.11 -0.18 1.17 -0.22 0.53 0.51
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.81
## beta[1] 1 0.52
## beta[2] 1 0.47
## beta[3] 1 0.39
## beta[4] 1 0.90
## gamma 1 0.73
## theta 1 0.76
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 4043 3612 8518 8605 3915 581 455
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 0.90 2.24 0.50 0.33 -0.26 -1.54 -0.84
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.32
## beta[1] 1801 0.21
## beta[2] 1 0.12
## beta[3] 1 0.22
## beta[4] 1 0.49
## gamma 1 0.08
## theta 1 0.34
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 3177 6390 8360 8768 3550 871 730
## Geweke diagnostics Z scores
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## -0.52 1.77 1.76 0.59 0.86 -0.65 -0.79
##
## Heidelberger and Welch's convergence diagnostic
## start pvalue
## beta[0] 1 0.13
## beta[1] 1 0.28
## beta[2] 1 0.40
## beta[3] 1 0.56
## beta[4] 1 0.13
## gamma 1 0.14
## theta 1 0.15
##
## Effective Sample Size
## beta[0] beta[1] beta[2] beta[3] beta[4] gamma theta
## 3544 5673 8501 9000 3650 724 730Table 5
cv<-function(gam,theta,mixing,ratio=FALSE)
{
if(mixing=="Gamma")
{
cv.star2=exp(lgamma(theta)+lgamma(theta-2/gam)-2*lgamma(theta-1/gam))-1
dcv.star2=exp(lgamma(theta)+lgamma(theta-2/gam)-2*lgamma(theta-1/gam))*(digamma(theta)+digamma(theta-2/gam)-2*digamma(theta-1/gam))
}
if(mixing=="InvGamma")
{
cv.star2=exp(lgamma(theta)+lgamma(theta+2/gam)-2*lgamma(theta+1/gam))-1
dcv.star2=exp(lgamma(theta)+lgamma(theta+2/gam)-2*lgamma(theta+1/gam))*(digamma(theta)+digamma(theta+2/gam)-2*digamma(theta+1/gam))
}
if(mixing=="InvGauss")
{
K1=besselK(x=1/theta,nu=-(1/gam+1/2)); K2=besselK(x=1/theta,nu=-(2/gam+1/2))
#Using the asymptotic expansion as theta->0
cv.star2.aux=sqrt(pi/2)*(exp(-1/theta)*theta^(1/2))*(K2/K1)*(1/K1)-1
aux1=sqrt(pi/2)*theta^(-3/2)*exp(-1/theta)/K1
aux2=(K2/K1)+besselK(x=1/theta,nu=-(2/gam-1/2))/K1-2*(K2/K1)*(1/K1)*besselK(x=1/theta,nu=-(1/gam-1/2))
dcv.star2.aux=aux1*aux2
cv.star2=ifelse(theta<0.0015,0,cv.star2.aux)
dcv.star2=ifelse(theta<0.0015,0,dcv.star2.aux)
}
if(mixing=="LogNormal")
{
cv.star2=exp(theta)-1
dcv.star2=exp(-2*log(gam)+theta)
}
aux0=(lgamma(1+2/gam)-2*lgamma(1+1/gam)); aux1=exp(aux0)-1
cv=sqrt(exp(aux0+log(cv.star2))+aux1)
if(ratio==FALSE){return(cv)}
if(ratio==TRUE){return(cv/sqrt(aux1))}
if(ratio=="WEI"){return(sqrt(aux1))}
}
cv.Weibull<-function(gam)
{
aux=sqrt(exp(lgamma(1+2/gam)-2*lgamma(1+1/gam))-1)
return(aux)
}
RcvGam = matrix(0, ncol = 3, nrow = 8)
for(i in 1:nrow(ParamTableGam))
{
aux = mcmc(cv(ChainGam[[i]]$gam,ChainGam[[i]]$theta,mixing="Gamma",ratio=TRUE))
RcvGam[i, 1] = median(aux)
RcvGam[i, 2:3] = HPDinterval(aux)
}
colnames(RcvGam) <- c("Ecv", "HPDEcvlower", "HPDEcvupper")
RcvIGam = matrix(0, ncol = 3, nrow = 8)
for(i in 1:nrow(ParamTableIGam))
{
aux = mcmc(cv(ChainIGam[[i]]$gam,ChainIGam[[i]]$theta,mixing="InvGamma",ratio=TRUE))
RcvIGam[i, 1] = median(aux)
RcvIGam[i, 2:3] = HPDinterval(aux)
}
colnames(RcvIGam) <- c("Ecv", "HPDEcvlower", "HPDEcvupper")
RcvIGauss = matrix(0, ncol = 3, nrow = 8)
for(i in 1:nrow(ParamTableIGauss))
{
aux = mcmc(cv(ChainIGauss[[i]]$gam,ChainIGauss[[i]]$theta,mixing="InvGauss",ratio=TRUE))
RcvIGauss[i, 1] = median(aux)
RcvIGauss[i, 2:3] = HPDinterval(aux)
}
colnames(RcvIGauss) <- c("Ecv", "HPDEcvlower", "HPDEcvupper")
RcvLN = matrix(0, ncol = 3, nrow = 8)
for(i in 1:nrow(ParamTableLN))
{
aux = mcmc(cv(ChainLN[[i]]$gam,ChainLN[[i]]$theta,mixing="LogNormal",ratio=TRUE))
RcvLN[i, 1] = median(aux)
RcvLN[i, 2:3] = HPDinterval(aux)
}
colnames(RcvLN) <- c("Ecv", "HPDEcvlower", "HPDEcvupper")
cat("Estimated ratios for each model \n")## Estimated ratios for each modelcbind(ParamTableGam, RcvGam)## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed Ecv HPDEcvlower
## 1 4 1.5 TruncExp 1 571629 2.391995 1.226446
## 2 1 1.5 TruncExp 1 492883 2.322744 1.190385
## 3 4 5.0 TruncExp 1 783694 6.986370 1.589244
## 4 1 5.0 TruncExp 1 540340 6.943675 1.448344
## 5 4 1.5 Pareto 1 88002 2.329099 1.087620
## 6 1 1.5 Pareto 1 420628 2.240850 1.040911
## 7 4 5.0 Pareto 1 976328 4.344738 1.146542
## 8 1 5.0 Pareto 1 330313 4.070118 1.159338
## HPDEcvupper
## 1 4.422618
## 2 4.235218
## 3 19.874946
## 4 19.789993
## 5 6.134094
## 6 5.767696
## 7 29.566035
## 8 27.916060cbind(ParamTableIGam, RcvIGam)## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed Ecv HPDEcvlower
## 1 4 1.5 TruncExp 1 215942 1.422416 1.243369
## 2 1 1.5 TruncExp 1 63169 1.407057 1.209227
## 3 4 5.0 TruncExp 1 521826 1.430164 1.237122
## 4 1 5.0 TruncExp 1 503616 1.403782 1.194593
## 5 4 1.5 Pareto 1 120486 1.399151 1.192560
## 6 1 1.5 Pareto 1 87275 1.364267 1.131292
## 7 4 5.0 Pareto 1 433560 1.422207 1.198182
## 8 1 5.0 Pareto 1 193080 1.400980 1.201202
## HPDEcvupper
## 1 1.554756
## 2 1.553597
## 3 1.556096
## 4 1.551413
## 5 1.551991
## 6 1.539696
## 7 1.570126
## 8 1.546437cbind(ParamTableIGauss, RcvIGauss)## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed Ecv HPDEcvlower
## 1 4 1.5 TruncExp 1 364452 1.671872 1.467143
## 2 1 1.5 TruncExp 1 772498 1.662514 1.452708
## 3 4 5.0 TruncExp 1 734878 1.676073 1.471130
## 4 1 5.0 TruncExp 1 972528 1.658172 1.435204
## 5 4 1.5 Pareto 1 740140 1.646229 1.395405
## 6 1 1.5 Pareto 1 200713 1.620578 1.352898
## 7 4 5.0 Pareto 1 377307 1.669073 1.436302
## 8 1 5.0 Pareto 1 247017 1.637123 1.390757
## HPDEcvupper
## 1 1.830358
## 2 1.831276
## 3 1.848358
## 4 1.838652
## 5 1.831929
## 6 1.817909
## 7 1.837317
## 8 1.829701cbind(ParamTableLN, RcvLN)## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed Ecv HPDEcvlower
## 1 4 1.5 TruncExp 1 97915 2.384842 1.533376
## 2 1 1.5 TruncExp 1 995633 2.238037 1.457315
## 3 4 5.0 TruncExp 1 320552 2.559325 1.689047
## 4 1 5.0 TruncExp 1 397726 2.431954 1.611196
## 5 4 1.5 Pareto 1 553224 2.283492 1.519334
## 6 1 1.5 Pareto 1 795749 2.075842 1.218768
## 7 4 5.0 Pareto 1 902245 2.434995 1.578776
## 8 1 5.0 Pareto 1 827927 2.247022 1.415542
## HPDEcvupper
## 1 3.169264
## 2 3.126312
## 3 3.388848
## 4 3.289663
## 5 3.159489
## 6 2.969993
## 7 3.291863
## 8 3.152258cat("Estimated true ratio w.r.t. Weibull model (Gamma model only with gam~Gamma(1,4))")## Estimated true ratio w.r.t. Weibull model (Gamma model only with gam~Gamma(1,4))aux = cv(ChainGam[[1]]$gam,ChainGam[[1]]$theta,mixing="Gamma") / cv.Weibull(ChainNone[[1]]$gam)
round(cbind(median(aux), HPDinterval(mcmc(aux))),2)## lower upper
## V1 2.08 1.08 3.74Model comparison
Pseudo ML
CDnone = lapply(ChainNone, FUN = RMWreg_CaseDeletion,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "None", BaseModel = "Weibull")
ParamTableNone$logPsML = sapply(CDnone, FUN = function(x) sum(x[,1]))
CDexp = lapply(ChainExp, FUN = RMWreg_CaseDeletion,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "Exponential", BaseModel = "Weibull")
ParamTableExp$logPsML = sapply(CDexp, FUN = function(x) sum(x[,1]))
CDgam = lapply(ChainGam, FUN = RMWreg_CaseDeletion,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "Gamma", BaseModel = "Weibull")
ParamTableGam$logPsML = sapply(CDgam, FUN = function(x) sum(x[,1]))
CDigam = lapply(ChainIGam, FUN = RMWreg_CaseDeletion,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "InvGamma", BaseModel = "Weibull")
ParamTableIGam$logPsML = sapply(CDigam, FUN = function(x) sum(x[,1]))
CDigauss = lapply(ChainIGauss, FUN = RMWreg_CaseDeletion,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "InvGauss", BaseModel = "Weibull")
ParamTableIGauss$logPsML = sapply(CDigauss, FUN = function(x) sum(x[,1]))
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
CDln = parLapply(cl, ChainLN, fun = RMWreg_CaseDeletion,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "LogNormal", BaseModel = "Weibull")
stopCluster(cl)
ParamTableLN$logPsML = sapply(CDln, FUN = function(x) sum(x[,1]))ML
RunRMWreg_logML <- function(Row, ParamTable, Chains,
Mixing, BaseModel,
Time, Event, DesignMat,
N, Thin, Burn, lambdaPeriod = 1)
{
require(RMWreg)
if(Mixing %in% c("None", "Exponential"))
{
logML = RMWreg_logML(Chains[[Row]], Time, Event, DesignMat,
Hyp1Gam = ParamTable$Hyp1Gam[Row], Hyp2Gam = ParamTable$Hyp2Gam[Row],
Mixing = Mixing, BaseModel = BaseModel,
N=N, Thin = Thin, Burn = Burn,
lambdaPeriod = lambdaPeriod)
}
else
{
logML = RMWreg_logML(Chains[[Row]], Time, Event, DesignMat,
PriorCV = as.character(ParamTable$PriorCV[Row]),
PriorMeanCV = ParamTable$PriorMeanCV[Row],
Hyp1Gam = ParamTable$Hyp1Gam[Row], Hyp2Gam = ParamTable$Hyp2Gam[Row],
Mixing = Mixing, BaseModel = BaseModel,
N=N, Thin = Thin, Burn = Burn,
lambdaPeriod = lambdaPeriod)
}
write.table(logML, paste0("logML_",Mixing,"_",Row,".txt"), row.names = FALSE, col.names = FALSE)
return(logML)
}
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableNone$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableNone)),
fun = RunRMWreg_logML,
ParamTable = ParamTableNone,
Chains = ChainNone,
Mixing = "None",
BaseModel = "Weibull",
Time = Time, Event = Event,
DesignMat = DesignMat,
N=N, Thin = Thin, Burn = Burn,
lambdaPeriod = 1))
stopCluster(cl)
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableExp$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableExp)),
fun = RunRMWreg_logML,
ParamTable = ParamTableExp,
Chains = ChainExp,
Mixing = "Exponential",
BaseModel = "Weibull",
Time = Time, Event = Event,
DesignMat = DesignMat,
N=N, Thin = Thin, Burn = Burn,
lambdaPeriod = 1))
stopCluster(cl)
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableGam$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableGam)),
fun = RunRMWreg_logML,
ParamTable = ParamTableGam,
Chains = ChainGam,
Mixing = "Gamma",
BaseModel = "Weibull",
Time = Time, Event = Event,
DesignMat = DesignMat,
N=N, Thin = Thin, Burn = Burn,
lambdaPeriod = 1))
stopCluster(cl)
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableIGam$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableIGam)),
fun = RunRMWreg_logML,
ParamTable = ParamTableIGam,
Chains = ChainIGam,
Mixing = "InvGamma", BaseModel = "Weibull",
Time = Time, Event = Event,
DesignMat = DesignMat,
N=N, Thin = Thin, Burn = Burn,
lambdaPeriod = 1))
stopCluster(cl)
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableIGauss$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableIGauss)),
fun = RunRMWreg_logML,
ParamTable = ParamTableIGauss,
Chains = ChainIGauss,
Mixing = "InvGauss", BaseModel = "Weibull",
Time = Time, Event = Event, DesignMat = DesignMat,
N=N, Thin = Thin, Burn = Burn,
lambdaPeriod = 1))
stopCluster(cl)
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableLN$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableLN)),
fun = RunRMWreg_logML,
ParamTable = ParamTableLN,
Chains = ChainLN,
Mixing = "LogNormal", BaseModel = "Weibull",
Time = Time, Event = Event, DesignMat = DesignMat,
N=N, Thin = Thin, Burn = Burn,
lambdaPeriod = 1))
stopCluster(cl)DIC
ParamTableNone$DIC <- sapply(ChainNone, FUN = RMWreg_DIC,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "None", BaseModel = "Weibull")
ParamTableExp$DIC <- sapply(ChainExp, FUN = RMWreg_DIC,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "Exponential", BaseModel = "Weibull")
ParamTableGam$DIC = sapply(ChainGam, FUN = RMWreg_DIC,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "Gamma", BaseModel = "Weibull")
ParamTableIGam$DIC = sapply(ChainIGam, FUN = RMWreg_DIC,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "InvGamma", BaseModel = "Weibull")
ParamTableIGauss$DIC = sapply(ChainIGauss, FUN = RMWreg_DIC,
Time = Time, Event = Event, DesignMat = DesignMat,
Mixing = "InvGauss", BaseModel = "Weibull")
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableLN$DIC = simplify2array(parLapply(cl, ChainLN, fun = RMWreg_DIC,
Time = Time, Event = Event,
DesignMat = DesignMat,
Mixing = "LogNormal",
BaseModel = "Weibull"))
stopCluster(cl)Summary
cat("\n RMW no mixing \n")##
## RMW no mixingParamTableNone## Hyp1Gam Hyp2Gam seed logPsML logML
## 1 4 1 307767 -1235.611 -1240.090
## 2 1 1 257673 -1235.667 -1238.788cat("\n RMW Exponential mixing \n")##
## RMW Exponential mixingParamTableExp## Hyp1Gam Hyp2Gam seed logPsML logML
## 1 4 1 372199 -1227.396 -1231.353
## 2 1 1 43825 -1227.506 -1230.724cat("\n RMW Gamma mixing \n")##
## RMW Gamma mixingParamTableGam## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed logPsML logML
## 1 4 1.5 TruncExp 1 571629 -1230.284 -1234.942
## 2 1 1.5 TruncExp 1 492883 -1230.294 -1233.972
## 3 4 5.0 TruncExp 1 783694 -1229.227 -1233.306
## 4 1 5.0 TruncExp 1 540340 -1229.285 -1232.612
## 5 4 1.5 Pareto 1 88002 -1230.321 -1235.412
## 6 1 1.5 Pareto 1 420628 -1230.487 -1234.411
## 7 4 5.0 Pareto 1 976328 -1229.537 -1234.149
## 8 1 5.0 Pareto 1 330313 -1229.655 -1233.316cat("\n RMW Inverse Gamma mixing \n")##
## RMW Inverse Gamma mixingParamTableIGam## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed logPsML logML
## 1 4 1.5 TruncExp 1 215942 -1230.997 -1231.651
## 2 1 1.5 TruncExp 1 63169 -1231.057 -1233.952
## 3 4 5.0 TruncExp 1 521826 -1230.985 -1235.228
## 4 1 5.0 TruncExp 1 503616 -1231.166 -1234.474
## 5 4 1.5 Pareto 1 120486 -1231.048 -1235.760
## 6 1 1.5 Pareto 1 87275 -1231.385 -1234.472
## 7 4 5.0 Pareto 1 433560 -1231.030 -1233.607
## 8 1 5.0 Pareto 1 193080 -1231.065 -1234.788cat("\n RMW Inverse Gaussian mixing \n")##
## RMW Inverse Gaussian mixingParamTableIGauss## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed logPsML logML
## 1 4 1.5 TruncExp 1 364452 -1228.718 -1233.280
## 2 1 1.5 TruncExp 1 772498 -1228.764 -1232.417
## 3 4 5.0 TruncExp 1 734878 -1228.781 -1233.141
## 4 1 5.0 TruncExp 1 972528 -1228.803 -1232.628
## 5 4 1.5 Pareto 1 740140 -1228.960 -1233.690
## 6 1 1.5 Pareto 1 200713 -1229.028 -1233.471
## 7 4 5.0 Pareto 1 377307 -1228.847 -1233.492
## 8 1 5.0 Pareto 1 247017 -1228.974 -1232.903cat("\n RMW Log-Normal mixing \n")##
## RMW Log-Normal mixingParamTableLN## Hyp1Gam PriorMeanCV PriorCV Hyp2Gam seed logPsML logML
## 1 4 1.5 TruncExp 1 97915 -1227.920 -1232.143
## 2 1 1.5 TruncExp 1 995633 -1228.143 -1231.776
## 3 4 5.0 TruncExp 1 320552 -1227.818 -1232.327
## 4 1 5.0 TruncExp 1 397726 -1228.007 -1232.309
## 5 4 1.5 Pareto 1 553224 -1228.010 -1232.904
## 6 1 1.5 Pareto 1 795749 -1228.567 -1232.309
## 7 4 5.0 Pareto 1 902245 -1227.877 -1232.262
## 8 1 5.0 Pareto 1 827927 -1228.286 -1231.900Figure 9
PosteriorSummary <- function(Column, ChainNone, ChainExp, ChainGam,
ChainIGam, ChainIGauss, ChainLN)
{
Median = c(unlist(lapply(ChainNone, function(x) median(x[,Column]))),
unlist(lapply(ChainExp, function(x) median(x[,Column]))),
unlist(lapply(ChainGam, function(x) median(x[,Column]))),
unlist(lapply(ChainIGam, function(x) median(x[,Column]))),
unlist(lapply(ChainIGauss, function(x) median(x[,Column]))),
unlist(lapply(ChainLN, function(x) median(x[,Column]))))
HPD = rbind(t(sapply(ChainNone, function(x) HPDinterval(x[,Column]))),
t(sapply(ChainExp, function(x) HPDinterval(x[,Column]))),
t(sapply(ChainGam, function(x) HPDinterval(x[,Column]))),
t(sapply(ChainIGam, function(x) HPDinterval(x[,Column]))),
t(sapply(ChainIGauss, function(x) HPDinterval(x[,Column]))),
t(sapply(ChainLN, function(x) HPDinterval(x[,Column]))))
return(cbind(Median, HPD))
}
plotHPD <- function(hpd, ...)
{
plotCI(x = hpd[,1], uiw = hpd[,3]-hpd[,1], liw = hpd[,1]-hpd[,2],
xaxt = "n", lwd = 2, xlab = "Mixing-Prior", gap = 1.5, ...)
}
param.name = c(expression(beta[0]), expression(beta[1]), expression(beta[2]), expression(beta[3]),
expression(beta[4]), expression(gamma), expression(theta))
titles = c(expression("Intercept"), expression("No impairments"),
expression("1 impairment"), expression("2 impairments"), expression("Birth weight"),"")
par(mfrow = c(6,1))
par(mar = c(7.5,11,6,1), mgp = c(6, 2, 0))
for(Column in 1:6)
{
PS = PosteriorSummary(Column, ChainNone_mcmc, ChainExp_mcmc, ChainGam_mcmc,
ChainIGam_mcmc, ChainIGauss_mcmc, ChainLN_mcmc)
if(Column == 1)
{
plotHPD(PS, main = titles[Column], ylab = param.name[Column],
cex.lab = 2.5, cex.axis = 2.5, cex.main = 4, ylim = c(2.5,4.8))
}
else
{
plotHPD(PS, main = titles[Column], ylab = param.name[Column],
cex.lab = 2.5, cex.axis = 2.5, cex.main = 4)
}
abline(v = c(2.5,4.5,8.5,12.5,16.5,20.5,24.5,28.5,32.5), lty = 1, lwd = 4)
abline(v = c(6.5,10.5,14.5,18.5,22.5,26.5,30.5,34.5), lty = 3)
axis(side = 1, at = c(1.5,3.5,6.5,10.5,14.5,18.5,22.5,26.5,30.5,34.5),
labels = c("None", "Exp", "Gamma-TExp", "Gam-Pare", "IGam-TExp", "IGam-Pareto",
"IGauss-TEx", "IGauss-Pare", "LN-TExp", "LN-Pareto"), cex.axis = 2)
if(Column == 1)
{
text(seq(5.5,36,by=2), rep(4.5,16),
rep(c(expression(paste("E(",c[v],")=1.5")),
expression(paste("E(",c[v],")=5"))),8), cex = 1.5)
}
}## Note: no visible binding for global variable 'v'
## Note: no visible binding for global variable 'v'
Figure 10
PlotModelComparison <- function(Row, ParamValue, Mains,
ParamTableNone, ParamTableExp, ParamTableGam,
ParamTableIGam, ParamTableIGauss, ParamTableLN, ...)
{
logBF = c(ParamTableNone$logML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableExp$logML[ParamTableExp$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" &
ParamTableGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableIGam$logML[ParamTableIGam$PriorCV == "TruncExp" &
ParamTableIGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableIGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "TruncExp" &
ParamTableIGauss$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableIGauss$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" &
ParamTableLN$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableLN$Hyp1Gam == ParamValue$Hyp1Gam[Row]]) -
ParamTableNone$logML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]]
logBF_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" &
ParamTableGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableIGam$logML[ParamTableIGam$PriorCV == "Pareto" &
ParamTableIGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableIGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "Pareto" &
ParamTableIGauss$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableIGauss$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" &
ParamTableLN$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableLN$Hyp1Gam == ParamValue$Hyp1Gam[Row]]) -
ParamTableNone$logML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]]
logPsBF = c(ParamTableNone$logPsML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableExp$logPsML[ParamTableExp$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" &
ParamTableGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "TruncExp" &
ParamTableIGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableIGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "TruncExp" &
ParamTableIGauss$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableIGauss$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" &
ParamTableLN$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableLN$Hyp1Gam == ParamValue$Hyp1Gam[Row]]) -
ParamTableNone$logPsML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]]
logPsBF_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" &
ParamTableGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "Pareto" &
ParamTableIGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableIGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "Pareto" &
ParamTableIGauss$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableIGauss$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" &
ParamTableLN$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
ParamTableLN$Hyp1Gam == ParamValue$Hyp1Gam[Row]]) -
ParamTableNone$logPsML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]]
plot(logBF, logPsBF, pch=c(4,8,0,5,1,2), cex=rep(3,6), bty = "n",
ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
ylim = c(0,9), xlim = c(0,9), main = Mains[Row], ...)
points(logBF_1, logPsBF_1, pch=c(15,18,16,17), cex=c(3,4.5,3,3), ...)
if(Row == 4)
{
legend('bottomright', c("Weibull", "RMWEXP", "RMWGAM", "RMWIGAM", "RMWIGAUSS", "RMWLN"),
pch = c(4,8,0,5,1,2), col = "black", cex = 1.5, bty = "n")
}
}
ParamValue = expand.grid("Hyp1Gam" = c(4, 1),
"PriorMeanCV" = c(1.5, 5))
Mains = rep(c(expression(paste(gamma,"~Gamma(",4,",",1,")")),
expression(paste(gamma,"~Gamma(",1,",",1,")"))),2)
par(mar=c(5,5,4,0.5))
par(mfrow=c(2,2))
for(Row in 1:4)
{
PlotModelComparison(Row, ParamValue, Mains,
ParamTableNone, ParamTableExp, ParamTableGam,
ParamTableIGam, ParamTableIGauss, ParamTableLN,
cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd = 2)
}
Figure 11
# Outlier detection
BFexp = RMWreg_BFoutlier(ChainExp[[1]], RefLambda = Event + (1-Event)*1/2,
Time, Event, DesignMat,
Mixing = "Exponential", BaseModel = "Weibull")
par(mar = c(5, 5, 4, 2) + 0.1)
plot(BFexp, type = "l", ylim = c(0, 2), main = "",
ylab = "Bayes Factors", xlab = "Patient",
cex.main = 3, cex.axis = 2, cex.lab = 2.5)
abline(h = 1, lty = 2)
Session information
sessionInfo()## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.3
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] compiler parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] KMsurv_0.1-5 data.table_1.9.6 gplots_3.0.1 coda_0.18-1
## [5] RMWreg_0.0.18 knitr_1.15.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.8 magrittr_1.5 lattice_0.20-34
## [4] stringr_1.1.0 highr_0.6 caTools_1.17.1
## [7] tools_3.3.2 grid_3.3.2 KernSmooth_2.23-15
## [10] htmltools_0.3.5 gtools_3.5.0 yaml_2.1.14
## [13] rprojroot_1.1 digest_0.6.10 bitops_1.0-6
## [16] evaluate_0.10 rmarkdown_1.2 gdata_2.17.0
## [19] stringi_1.1.2 backports_1.0.4 chron_2.3-47